Non linear transport in drift-diffusion equations under magnetic field 
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We analyze numerically and analytically the non linear transport properties of a drift-diffusion 
equation in presence of a magnetic field and of a disorder potential. For a wide range of parameters 
this model exhibits a plateau where the drift velocity is almost independent on the applied electric 
field. This behavior has strong similarities with the zero differential resistance states observed 
experimentally in high mobility two dimensional systems. Performed numerical simulations are in a 
good global agreement with the developed analytical theory even if the later leads to overestimated 
negative differential resistance values. 
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The study of out of equihbrium transport in high pu- 
rity two dimensional electron systems has lead to the 
discovery of several surprising effects t lpSt and stimu- 
lated an intense theoretical research [9l4l6|. Probably 
the most striking phenomenon observed so far is the for- 
mation of zero-resistances states under microwave irradi- 
ation at certain ratios between microwave and cyclotron 
frequencies 17H23l| . It is widely argued that in a zero re- 
sistance state the electron gas undergoes a transition to 
a state where the Hall current splits into several domains 
IS 16 1 ■ Zero resistance is then explained as the 
result of a cancellation between the contributions of the 
different domains. Although a microscopic theory de- 
scribing the domain formation remains unavailable, the 
domain concept was used to interpret many other ex- 
periments. Thus when zero differential resistance states 
(ZDRS) were first reported, it was assumed that they 
were also a manifestation of this domain physics [24L 
[2^. However recent experiments on ZRS in liquid He- 
lium show that the steady state of the surface electrons 
in this regime is not a domain state but an accumulation 
of the electrons on the edges of the electron cloud Jl] . 
This and recent alternative theoretical proposals which 

stimulated 



do not rely on the domain picture |32h34 
us to see if an alternative explanation was possible for 
ZDRS as well. Thus in this letter we propose a simple 
transport model which produces an effect very similar to 
ZDRS without appealing to a domain formation picture. 



We will start by introducing the transport model em- 
ployed in our investigations. A carrier in a two dimen- 
sional electron gas moves under the influence of a disorder 
potential U{r) — / ^|^J7(q)e*'''', electron-electron and 
electron-phonon interactions. According to results from 
classical dynamics [s^, an electron will follow adiabati- 
cally all the fluctuations of the potential on a wavelength 
2n/q much larger than the Larmor radius when placed 
under a magnetic field. Thus we separate the disorder po- 
tential into two independent contribution, the long wave- 
length part which we henceforth call U (r) for simplic- 
ity and a short wavelength contribution. Our next as- 
sumption is that the short- wavelength disorder, electron- 



phonon and electron-electron interactions give rise to a 
"microscopic" longitudinal mobility and to a "micro- 
scopic" diffusion constant D. Under these assumptions 
the simplest model to describe electron transport is a 
drift-diffusion equation: 



dtP = div (Aferad U - E]P) + DAP 



(1) 



where P(r) is the probability distribution of the elec- 
tron, /t is the mobility tensor with diagonal and off- 
diagonal components ^xx and ±iixy respectively. We 
have also introduced the applied static electric field: 
E. To specify the model completely we need to de- 
fine the correlator A(r — r') =< J7(r)[/(r') >. In two 
dimensions a possible form for its Fourier transform is 
A{q) = exp 2Q^) ' '^here A gives the disorder am- 
plitude and the scale determines the typical wave- 
length for the fluctuations of the potential. Although A 
is not expected to depend on the magnetic field, the scale 
and the microscopic mobility fixx depend the Lar- 
mor radius since it gives the characteristic length-scale 
where the adiabatic theory starts to apply. In the follow- 
ing we do not try to derive quantitative predictions for 
these parameters, instead we focus on the general prop- 
erties of Eq. ([T]) and show that it produces a non-linear 
transport behavior very similar to ZDRS. 

Drift diffusion equations in random media has been 
studied in the past in several contexts including fluid dy- 
namics and reaction kinetics in cells 3B"40[. The pres- 
ence of a random potential U{r) is expected to change 
the microscopic diffusion rate and mobility to effective 
large scale values D* and /i*^. The originality of Eq. ([T]) 
comes from the non diagonal mobility tensor fl (in high 
mobility samples fj,xy 3> fJ,xx at even moderate magnetic 
fields i? ~ 0.1 Tesla) and the presence of an applied elec- 
tric field E which can distort the disorder potential U (r). 
As a result we look not only for the effective mobility /i*^ 
but for the full dependence of the longitudinal drift ve- 
locity Vx on the the applied electric field Ex- ZDRS cor- 
responds to a regime where the drift velocity Vx does not 
depend on the bias Ex- Indeed in a Hall bar the current 
density perpendicular to the channel vanishes far away 
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FIG. 1: Dependence of the drift velocity as a function of 
the appUed field Ex for several values of the microscopic dif- 
fusion coefficient D (in dimensionless units). The dynamics 
of the probe trajectories was followed for time 10^ Hxy^Q^ 
and Nq = 100 disorder Harmonics were used in the real- 
ization of the disorder potential. The inset show the equi- 
potential curves for a potential with A'^, = 3 under electric 
fields Ex/iQX) = 0.25, 0.5 (left/right). The thick red seg- 
ments indicate the thickness of the bundle of delocalized equi- 
potentials (approximately proportional to Ex). 



from the current injecting electrodes, which leads to an 
expression for the longitudinal field = Vx{E±_) / as 
a function of the Hall field E±^ = „ where / the in- 

Kh Wc 

jected current, Rh is the Hall resistance and Wc is the 
channel width (simulations on the basis of Eq. ([T]) and 
experiments indicate that the Hall mobility fx^y and Hall 
resistance Rh remain almost unchanged by the applied 
electric field). 

The average drift velocity for the model Eq. ([T]) can 
be determined by integrating the Langevin equations of 
motion v = /i (E — grad U)+£, for several probe particles. 
In this equation v is the particle velocity and ^ is a noise 
term whose amplitude is chosen to create the microscopic 
diffusion rate D. In order to simulate the random poten- 
tial, we have taken U{x, y) = ' A]^ cos(kr -I- 0k) 
where the amplitudes and wavevectors k were taken 
from a Gaussian distribution with standard deviations 
A-\/2 and Q respectively (the phases (/)k are uniformly 
distributed in (0,27r)) . The simulations were performed 
using a Runge-Kutta integration method and typically 
Nk = 100 harmonics were used for the realization of 
the disorder potential. On Fig. 1 we represent the de- 
pendence of Vx on the applied field, for several values 
of the microscopic diffusion rate D. In absence of the 
disorder potential the drift velocity would be given by 
Vxx — fJ-xxEx, this relation is fairly accurate for high val- 
ues of D/{iJ,xyX). However when the diffusion rate D is 
lowered, the dependence Vx{Ex) becomes highly nonlin- 
ear and starts to exhibit a plateau where Vx depends very 
weakly on the applied field Ex ■ For some values of D our 



simulations suggest a weak decrease of Vxx{Ex) (around 
5%) in this region however the obtained behavior is still 
remarkably flat. We also noticed than the small drop in 
drift velocity became more pronounced for a small num- 
ber of Harmonics, thus the plateau may become com- 
pletely flat in the limit Nk oo. 

The origin for this plateau can be understood by rep- 
resenting the equipotentials in presence of a moderate 
applied electric field Ex in the x direction. The po- 
tential landscape is formed by closed equipotential is- 
lands which are separated by narrow streams of delocal- 
ized equipotentials which flow perpendicularly (on av- 
erage) to the applied field. The thickness W of these 
streams grows linearly with the applied field Ex and can 
be estimated as ~ ExQ~'^\~^ Thus the time 

Ts ~ W / {pixxEx) ~ Q~^X~^ / fjLxx to drift across the 
streams is independent of the applied field. After crossing 
such a stream a carrier is transferred to another island. In 
the process it has drifted in the direction of the applied 
field by a quantity Q^^, leading to a saturation value 
for the drift velocity: Vsat ^ Q^^ /Tg ^ fJ-xxXQ^^. This 
reasoning is plausible as long as the equipotentials are 
weakly distorted by the applied field Ex <^ XQ; for higher 
values of Ex we expect to recover Vx = fJ-xxEx as in the 
absence of a disorder potential. Although the above ar- 
gument yields the correct scaling (as compared to Fig. 1), 
the derivation is arguable in many ways. Thus to pursue 
analytical investigations we have employed more stan- 
dard field theoretical and renormalization group tech- 
niques, following closely the presentation from Ref. [39| . 

The drift velocity can be obtained from the asymp- 
totic properties of the disorder averaged Green func- 
tion for Eq. dl]) in the hmit k — > 0, indeed G(k)^^ = 
i{kxVx + kyVy) -\- 0{k^). In the Green-function repre- 
sentation, the perturbation theory series in the power 
of disorder potential A can be conveniently represented 
through Feyman diagrams. The resulting diagrammatic 
expansion rules are as follow: full lines represent the bare 
Green function in absence of the disorder potential 

dotted lines carry a factor A(q), and the interac- 
tion vertex is (k, fiq) where k and q are the incom- 
ing electron/" photon" wavenumbers respectively. The 
wavenumbers are conserved at each vertex and wavenum- 
ber is integrated around closed loops with a factor 
d2g/(2^)2. 

These diagrammatic rules lead us to the following one- 
loop expression for the free energy (the corresponding 
diagram is shown on Fig. 2) 

^A(g)(k,/iq)(k + q, Aq)Go(k + q) 

(3) 
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we remind that the free energy E(k) is connected to the 
Green function through the Dyson equation G~^(k) = 
Gq^{\s.) — S(k). The correction to the velocity is 
obtained by keeping only the linear terms in from 
Eq. ^ and neglecting all second order terms in the quan- 
tity Hxx/^J'xy ^ 1- Under this approximation the in- 
tegral over momentum q can be computed analytically 
and leads to: 



Vx — fJ'XxEx + A 



2 f^xxQ j^ff^xy-^x^ 



Ex 



-rC- 



DQ 



(4) 



where the function R(x) is defined as: 



R{x) = 2- ^e"'/4 [x^KoixyA) + (2 - x^)Ki{xy4)] 

with Kq and Ki designating the modified Bessel func- 
tions. Although this expression compares favorably with 
the numerical simulations in the limit Hxy^/D — 0, it 
introduces a pronounced negative differential resistance 
state at strong A which is not present in the numerical 
simulations (see Fig. 2). 

The failure of one loop perturbation theory suggests 
that higher order nonlinear terms in the disorder poten- 
tial amplitude must also be taken into account. A pos- 
sible approach is to use the self-consistent Born approx- 
imation where the bare Green function Go (k) in Eq. ^ 
is replaced by the full propagator G(k). In order to solve 
the self-consistent equation, we have searched for a solu- 
tion of the form G(k)"^ = i{vxkx + Vyky) - D^k^ - D*ky 
and kept only first and second order terms in wavenum- 
ber that stem from Eq. Then we determined £>* and 
D* numerically using the approximation Vy ~ ^fJ-xyEx 
and Vx ^ Vy and used these value to find Vx- Under this 
approximation, we find effective diffusion rates D* and 
D* much stronger than the microscopic value D. This is 
a consequence of the enhancement of the diffusion rate 
on large length-scales due to the presence of the disorder 
potential. However for small D the underlying micro- 
scopic dynamics is almost ballistic at small length-scales, 
and this fact is missed when the same effective diffusion 
constants D* and D* is used for all length-scales. As 
a result the self-consistent Born approximation leads to 
a smooth weakly non-linear increase similar to the lin- 
ear response theory for high D (see Fig. 2) and does not 
reproduce the plateau in Vx{Ex)- 

The failure of the self-consistent Born approximation, 
suggests to use a renormalization group (RG) method 
where D and fixx are allowed to vary as a function of the 
length scale A~^. In the absence of magnetic and electric 
fields, [11] showed that RG gives an accurate effective 
diffusion rate. Thus we adapt the calculations developed 
therein to the case of magnetotransport in the regime 
fJ-xy S> fJ-xx and in the limit of vanishing electric fields Ex ■ 
In this case the functional form of the Green function at 
a length scale A~^ is simply GA(k)^^ — ~D{A)k^. As in 
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FIG. 2: Top figures: a) diagram for the one loop expression 
of the self energy leading to Eq. and b) represents the 
diagram for the vertex correction. The bottom curves show 
a comparison between the numerical calculation of Vx{Ex) 
and the different analytical approximations described in the 
text (one loop perturbation theory, self consistent Born ap- 
proximation and renormalization group theory) for the case 
of D/{ii^y\) = 0.05 in Fig. 1. 



perturbation theory, the one-loop diagram allows to com- 
pute the change in self-energy Sa when the momentum 
scale is reduced from A to A — (5 A: 



EA_A-A(k) - EA(k) = 



2_J^;xy_,2 



A(A) 
47r 



MK 



(5) 



The Dyson equation gives a relation between the change 
in the self-energy and the change in the diffusion con- 
stant -dT,/^/dA = k'^dD{A)/dA. We are thus lead to 
a first order differential equation on D{A), which after 
integration gives: 



2 / ' 



D* = D(0) (6) 



This expression for D* is in good agreement with the nu- 
merical simulations in a wide range of parameters. The 
regime D ^ deserves special attention; in this limit 
electrons are likely to become trapped around minima of 
the disorder potential, which results in vanishing trans- 
port coefficients D* , fi*^, fi*y. Since the RG calculation 
did not take this possibility into account, its domain of 
validity is restricted by 13 3> Xfixx where trapping is not 
expected. 

The effective mobility /i*^. can be deduced by includ- 
ing the vertex correction represented on Fig. 2 in the 
RG procedure. Indeed in the linear response regime, 
fixx appears in the diagrams only through the inter- 
action vertex. The ffow equation for i^xxiA) is then 
d^ixx{A) = fixx{A)D{A)-'^dD{A), which leads fj.*^^ = 
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fJ-xx{0) = D*fixx/D. This result can also be understood 
intuitively from the Einstein-relation, indeed the ratio 
D{A) I /i^j. (A) determines an effective carrier temperature 
which is not expected to depend on the length scale A~^. 
This expression for the mobility combined with Eq. ([5]) 
leads to an expression for the mobility which compares 
favorably with our numerical simulations in the regime 
of weak electric fields. 

In order to generalize the RG calculation to the case 
of finite electric fields several obstacles arise. The drift 
velocity is determined by the asymptotic properties of 
the Green function and not only by the renormalized in- 
teraction vertex; moreover the bare propagator becomes 
anisotropic which makes separation of momentum space 
in isotropic shells A — J A < r < A inappropriate. In 
front of these difhculties we have adopted the following 
simplifying anzats, we assumed that the renormalization 
flow equation for the drift velocity (A) was of the form 

v^{K) dA D{A) dA "\D{A)a) ^' 

The proposed form is a generalization of the flow equa- 
tion for the mobility Hxx{A). It incorporates the de- 
pendence on the electric field through the scaling func- 
tion Sv which depends on the dimensionless parameter 
IJ-xyEx/ {D{A)A). The choice of this quantity is governed 
by dimensional analysis and by similarity to Eq. (j4]) . The 
connection to Eq. ^ can be exploited even further, by 
requiring that Eq. ([7]) reduces to Eq. (|4]) in the limit of 
weak disorder A — > 0. This allows to find the functional 
form of Sy{x): 

The drift velocity Vx can now be obtained through nu- 
merical integration of Eq. (O using the initial condition 
Vx{A = oo) = fixxEx and Eqs. (I6]),([8|). The obtained val- 
ues for the drift velocity are in a good agreement with 
numerical simulations (see Fig. 2). They reproduce with 
a good accuracy the behavior in the limit of low and high 
electric fields, however the above approximation also pro- 
duces a clear negative differential resistance which is not 
present in the simulations. 

Several factors may explain the observed discrepancy. 
The change of the diffusion rate by the applied electric 
field is not captured in Eq. (|6]). Under bias the diffu- 
sion also becomes anisotropic which requires to intro- 
duce two diffusion coefficients Dx{A) and Dy{A) as in 
the self-consistent Born approximation. Moreover the 
change in Vx (A) deduced from the one loop expression of 
the self-energy and the change in ^xx (A) from the vertex- 
correction must be treated separately. If the isotropic RG 
procedure is still valid in presence of anisotropy, it may 
be possible to derive coupled equations on the flow of the 
four quantities: Dx{A), Dy{A), ^xx{A) and Vx{A), how- 
ever the results obtained from the anzats Eq. ([7]) already 



show a sufficient agreement with the simulations to give 
confidence in our numerical data and in the truthfulness 
of the observed the plateau in Vx{Ex). More advanced 
RG calculations will not necessarily bring a physical un- 
derstanding of the plateau origin, and other analytical 
methods could be more successful for this purpose [i^. 
Thus we did not try to push the accuracy of the RG 
calculations further. 

So far we did not discuss the role played by the for- 
mation of Landau levels under magnetic field. The first 
experimental reports [2J] indicated that ZDRS appeared 
only at magnetic fields corresponding to the maximum of 
a Shubnikov-de Haas (SdH) oscillation, while the devel- 
oped theory seem to suggest that they could appear at 
any magnetic field. We interpret this apparent contra- 
diction with the experiments as follow; the experiments 
were conducted in a magnetic field regime where the lon- 
gitudinal resistance vanished at the minimum of the SdH 
oscillations. It is known that in this regime the drop 
of the Hall voltage along the channel maybe strongly 
nonuniform outside of the SdH maxima [43l . 13] which 
could lead to an effective smoothing of the ZDRS. We 
note that the recent observation of ZDRS at high filling 
factors [4^ before the onset of SdH oscillations seems to 
support this interpretation. 

To summarize we have shown through numerical sim- 
ulations that model Eq. ([T|) leads to a non linear trans- 
port behavior very similar to a zero-differential-resistance 
state without invoking an instability argument. Even 
if our simulations are consistent with a weak residual 
negative-differential resistance, this trend is strongly ex- 
aggerated in some of the analytical approximations which 
we employ, specially if only the first order disorder con- 
tribution to the self energy is kept. Approximate renor- 
malization group calculations were the most successful 
to reach an agreement between simulation and analyt- 
ical theory, however even in this case the plateau was 
not reproduced correctly. For a quantitative comparison 
with experiments the parameters of the model must be 
estimated from microscopic transport theory which will 
be the focus of further attention. Moreover edge effects 
are probably also important, since it is likely that edge- 
trajectories are stabilized against the fiuctuations of the 
disorder potential by a static electric field perpendicu- 
lar to the edge. To conclude it is interesting that the 
proposed simple model already yields a rich non-linear 
physics and can potentially rise an interesting challenge 
for methods of field theory and non-linear science. 

We acknowledge M. Schindler and D.L. Shepelyansky 
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terest in this work and RIKEN for hospitality during the 
writing of the manuscript. 
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